library(dplyr)
library(INLA)
library(inlabru)
library(sf)
library(terra)
# load some libraries to generate nice map plots
library(scico)
library(ggplot2)
library(patchwork)
library(mapview)
library(tidyterra)Libraries to load:
The data
In this practical, we will revisit the data on the Pacific Cod (Gadus macrocephalus) from a trawl survey in Queen Charlotte Sound. The pcod dataset is available from the sdmTMB package and contains the presence/absence records of the Pacific Cod during each surveys along with the biomass density of Pacific cod in the area swept (kg/Km\(^2\)). The qcs_grid data contain the depth values stored as \(2\times 2\) km grid for Queen Charlotte Sound.
The dataset contains presence/absence data from 2003 to 2017. Lets consider year 2003 for now.
library(sdmTMB)
pcod_df = sdmTMB::pcod
qcs_grid = sdmTMB::qcs_gridThen, we create a sf object and assign the right coordinate reference to it:
pcod_sf = st_as_sf(pcod_df, coords = c("lon","lat"), crs = 4326)
pcod_sf = st_transform(pcod_sf,
crs = "+proj=utm +zone=9 +datum=WGS84 +no_defs +type=crs +units=km" )We convert the covariates into a raster and assign the same coordinate reference:
depth_r <- rast(qcs_grid, type = "xyz")
crs(depth_r) <- crs(pcod_sf)
Spatio-temporal LGOCV comparison
Model fitting
Now lets compare two different space-time models using LGOCV and some information criteria metrics. The general model structure is given by:
\[ \begin{aligned} y(s,t)|\eta(s,t)&\sim\text{Binom}(1, p(s,t))\\ \eta(s,t) &= \text{logit}(p(s,t)) \\ \end{aligned} \]
- Model 1 - time iid effect We consider a separable space-time model with a linear predictors given by:
\[ \eta(s,t) = \beta_0 + f_1(\text{depth}(s)) + f_2(t) + \omega(s) \]
\(f_1(\text{depth}(s))\) is a smooth covariate effect of depth (modeled using a RW2 model)
\(f_2(t)\) is an IID effect of time
\(\omega(s)\) is Matérn random field.
The first step is to define the mesh and the spde model
Construct the mesh and the SPDE model
mesh = fm_mesh_2d(loc = pcod_sf,
cutoff = 1,
max.edge = c(10,20),
offset = c(5,50),
crs = st_crs(pcod_df))
spde_model = inla.spde2.pcmatern(mesh,
prior.sigma = c(1, 0.5),
prior.range = c(100, 0.5))
create time index and the grouped variable
To use the RW2 model the covariate has to be groupes:
depth_r$depth_group = inla.group(values(depth_r$depth_scaled))we also define a time index from 1 to in the data frame
pcod_sf = pcod_sf %>%
mutate(time_idx = match(year, c(2003, 2004, 2005, 2007, 2009, 2011, 2013, 2015, 2017)),
id = 1:nrow(.)) # Observation id for CV##Implement Model 1##
Note
Note that there are some survey locations in certain years that fall outside the depth raster region. inlabru will input these missing covariate values using the nearest available value. This can be computationally expensive, but you can avoid it by supplying a raster layer that encompasses all of your data points (e.g., by pre-imputing these missing values with your preferred method of choice).
One way of doing this in the inlabru framework is to use the bru_fill_missing() function:
# Select the raster of interest
depth_orig = depth_r$depth_group
re <- extend(depth_orig, ext(depth_orig)*1.05)
# Convert to an sf spatial object
re_df <- re %>% stars::st_as_stars() %>% st_as_sf(na.rm=F)
# fill in missing values using the original raster
re_df$depth_group = bru_fill_missing(depth_orig,re_df,re_df$depth_group)
# rasterize
depth_filled <- stars::st_rasterize(re_df) %>% rast()
plot(depth_filled)
We have now fit the model and want to check the results.
- Model 2 - spatiotemporal field We consider a separable space time model with a linear predictor given by:
\[ \eta(s,t) = \beta_0 + f_1(\text{depth}(s)) + \omega(s,t) \]
- \(f_1(\text{depth}(s))\) is a smooth covariate effect of depth (RW2)
- \(\omega(s,t)\) is a space-time Matérn spatial field with AR1 time component \[ \omega(s,t) = \phi\ \omega(s,t-1) + \epsilon(s),\qquad \epsilon(s)\sim\text{GF}(\sigma_{\epsilon},\rho_{\epsilon}) \]